Trajectory Analysis
Thyme Lab
2023-01-09
Trajectory analysis of zebrafish pallium clusters using Monocle3.
1 Setup
#install version 4.1.1 of seurat
#remotes::install_version("Seurat", version = "4.1.1")
#install seurat wrappers
#devtools::install_github("satijalab/seurat-wrappers")
#BiocManager::install(c("Biobase", "SingleCellExperiment", "batchelor", #"BiocGenerics", "DelayedArray", "DelayedMatrixStats", "limma", "S4Vectors", #"SummarizedExperiment", "pcaMethods"))
#install dynverse
#devtools::install_github("dynverse/dyno")
#load packages
library(tidyverse) #tidyverse dataframe manipulation
library(Seurat) #scRNAseq analysis
library(monocle3) #pseudotype analysis
library(SeuratWrappers) #convert Seurat object to monocle3
library(patchwork) #graph layout
set.seed(1234)2 Subset Seurat Object
2.1 Load and subset Seurat Object
Subset original forebrain.integrated Seurat object to include Comm. Pa Precursors, Pallium IN01, Pallium IN02, Pallium 02, Pallium 03, Pallium 05, Pallium 06, and Pallium 07.
#load integrated object
load(file = "forebrain_all_final.RObj")
#rename two clusters
forebrain.integrated <- RenameIdents(object = forebrain.integrated, `UnK` = "NK-like")
forebrain.integrated <- RenameIdents(object = forebrain.integrated, `Astrocytes` = "Astrocyte-like")
#rename clusterinterp2
forebrain.integrated$clusterInterp2 <- recode_factor(forebrain.integrated$clusterInterp2, "Astrocytes" = "Astrocyte-like", "Pallium_09" = "Pallium 09", "PoA_02" = "PoA 02", "PoA_01" = "PoA 01", "Subpallium_03" = "Subpallium 03", "Pallium_IN03" = "Pallium IN03", "Progenitor_01" = "Progenitor 01", "Progenitor_02" = "Progenitor 02", "UnK" = "Nk-like", "Cycling progenitors_02" = "Cycling progenitors 02", "Subpallium_IN" = "Subpallium IN", "Olfactory Bulb" = "Olf. Bulb", "Cycling progenitors_01" = "Cycling progenitors 01", "Committed_pallium_precursors" = "Comm. Pa Precursors", "Pallium_IN01" = "Pallium IN01", "Pallium_IN02" = "Pallium IN02", "Pallium_01" = "Pallium 01", "Pallium_02" = "Pallium 02", "Pallium_03" = "Pallium 03", "Pallium_04" = "Pallium 04", "Pallium_05" = "Pallium 05", "Pallium_06" = "Pallium 06", "Pallium_07" = "Pallium 07", "Pallium_08" = "Pallium 08", "Subpallium_01" = "Subpallium 01", "Subpallium_02" = "Subpallium 02", "Subpallium_03" = "Subpallium 03", "Subpallium_04" = "Subpallium 04", "Subpallium_05" = "Subpallium 05", "Subpallium_06" = "Subpallium 06", "Subpallium_07" = "Subpallium 07", "Subpallium_08" = "Subpallium 08", "Committed_subpallium_precursors" = "Comm. SPa Precursors")
levels <- c('dHabenula', 'vHabenula', 'Microglia', 'OPC', 'Oligodendrocytes', 'Olfactory Bulb', 'Progenitor_02', 'Progenitor_01','NK-like', 'Epithelial/Endothelial Cells', 'Cycling progenitors_01', 'Cycling progenitors_02', 'Committed_pallium_precursors', 'Committed_subpallium_precursors', 'Pallium_IN01', 'Pallium_IN02', 'Subpallium_IN', 'Subpallium_03', 'Pallium_01', 'Pallium_02', 'Pallium_03', 'Pallium_04', 'Pallium_05', 'Pallium_06', 'Pallium_07', 'Pallium_08', 'Pallium_09', 'Subpallium_01', 'Subpallium_02', 'PoA_02', 'Subpallium_04', 'Subpallium_05', 'Subpallium_06', 'Subpallium_07', 'Subpallium_08', 'PoA_01', 'Astrocyte-like','Pallium_IN03')
Idents(forebrain.integrated) <- factor(x = Idents(forebrain.integrated), levels = levels)
#clusters that we want to use in the trajectory
clusterstokeep <- c("Comm. Pa Precursors", "Pallium IN01", 'Pallium 05', 'Pallium 06', 'Pallium 07', "Pallium IN02", "Pallium 02", "Pallium 03")
#subset integrated object by cluster
forebrain.integrated <- subset(forebrain.integrated, subset= clusterInterp2 %in% clusterstokeep)
#colors for selected clusters
cols <- c('lightsteelblue2', 'lightcyan1', 'darkslategray1', 'blue1', 'darkturquoise','mediumspringgreen','lightskyblue','royalblue1')
cols <- lapply(cols, adjustcolor, alpha= 1.0)
#UMAP plot
DimPlot(forebrain.integrated, cols = cols, label = F, pt.size = 0.6)
#make a copy to use later
forebrain.integrated.copy <- forebrain.integrated3 Trajectory anlaysis using Monocle
3.1 Convert Seurat object to CDS & cluster
#change default assay to integrated because it has UMAP data
DefaultAssay(forebrain.integrated) <- "integrated"
#convert seurat object to CDS
cds <- as.cell_data_set(forebrain.integrated)
#cluster and partition CDS
cds <- cluster_cells(cds)
#plot clusters and partitions
p1 <- plot_cells(cds, color_cells_by = "cluster", show_trajectory_graph = FALSE)
p2 <- plot_cells(cds, color_cells_by = "partition", show_trajectory_graph = FALSE)
wrap_plots(p1, p2)
Clustering looks similar to clustering conducted by Seurat. All cells
are in partition 1.
3.2 Learn graph & plot trajectory
#converting to seurat object and back prevents future error
forebrain.integrated <- subset(as.Seurat(cds, assay = NULL))
cds <- as.cell_data_set(forebrain.integrated)
#fit a principal graph within each partition
cds <- learn_graph(cds, use_partition = FALSE, verbose = FALSE)##
|
| | 0%
|
|======================================================================| 100%#plot trajectory
plot_cells(cds,
color_cells_by = "cluster",
label_groups_by_cluster=FALSE,
label_leaves=TRUE,
label_branch_points=TRUE)
3.3 Choose root & plot based on pseudotime
Root is in the part of the committed pallium precursors cluster that is closest to cycling progenitors.
#choose root
cds <- order_cells(cds, reduction_method = "UMAP", root_cells = colnames(cds[, clusters(cds) == 9]))
#plot based on pseudotime
plot_cells(cds,
color_cells_by = "pseudotime",
group_cells_by = "cluster",
label_cell_groups = FALSE,
label_groups_by_cluster=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
label_roots = FALSE,
cell_size = 0.6)
3.4 Add expression data to CDS
The integrated Seurat doesn’t have gene expression data for all genes. This step makes a new CDS object from the RNA assay of forebrain.integrated, then adds the expression data to the CDS object. The new CDS object should have trajectory constructed using the 2000 variable features but expression data from the RNA assay.
#change default assay to RNA
DefaultAssay(forebrain.integrated.copy) <- "RNA"
#make new monocle object from RNA assay with expression data for all genes
RNA_from_seurat <- as.cell_data_set(forebrain.integrated.copy)
#make a copy in case I mess up CDS
cds_test <- cds
#add expression data to monocle object
cds_test@assays <- RNA_from_seurat@assays
cds_test@elementMetadata <- RNA_from_seurat@elementMetadata
cds_test@int_elementMetadata <- RNA_from_seurat@int_elementMetadata
cds_test@rowRanges <- RNA_from_seurat@rowRanges3.5 Find genes that vary over pseudotime
#find genes that vary over pseudotime
#don't knit because it takes a long time
#cds_graph_test_results <- graph_test(cds_test,
# neighbor_graph = "principal_graph",
# cores = 8)
#export csv with genes (import this later to avoid running graph_test)
#write.csv(cds_graph_test_results, "cds_graph_test_results.csv")3.6 Graph genes that vary over pseudotime
#import genes from differential expression analysis - or use code above
cds_graph_test_results <- read.csv("cds_graph_test_results.csv", header=TRUE)
#makes genes into rownames of cds_graph_test_results
rownames(cds_graph_test_results) <- cds_graph_test_results$X
#makes a list of the genes with q < 0.05
deg_ids <- rownames(subset(cds_graph_test_results[order(cds_graph_test_results$morans_I, decreasing = TRUE),], q_value < 0.05))
#add gene names to CDS and CDS_test
rowData(cds)$gene_short_name <- row.names(rowData(cds))
rowData(cds_test)$gene_short_name <- row.names(rowData(cds_test))
#import genes available in AGETAZ database
agetazgenes <- read.csv("genesinagetaz.csv", header=FALSE, stringsAsFactors=FALSE)
agetazgenes <- as.character(agetazgenes[,1])
#genes to look up in AGETAZ
cds_graph_test_results$genenames <- rownames(cds_graph_test_results)
deg_ids_insitu <- cds_graph_test_results %>% filter(genenames %in% intersect(agetazgenes, deg_ids))
deg_ids_insitu <- rownames(subset(deg_ids_insitu[order(deg_ids_insitu$morans_I, decreasing = TRUE),], q_value < 0.05))
#genes in AGETAZ database that vary over pseudotime
deg_ids_insitu## [1] "her15.2" "otpb" "otpa"
## [4] "her4.1" "eomesa" "ybx1"
## [7] "her4.4" "hmga1a" "neurod4"
## [10] "ets2" "bhlhe23" "jun"
## [13] "snap25a" "her2" "bhlhe22"
## [16] "zic3" "nr4a1" "egr1"
## [19] "neurod6a" "insm1a" "hmgb2b"
## [22] "emx3" "nr2f2" "irx1a"
## [25] "nr2f1a" "fosb" "tbr1b"
## [28] "pbx3b" "egr2b" "tcf7l2"
## [31] "lhx5" "her8.2" "lhx9"
## [34] "hmgb2a" "barhl2" "nhlh2"
## [37] "scrt1b" "tshz1" "nfia"
## [40] "neurog1" "zic4" "id2a"
## [43] "fosl2" "phactr3b" "si:ch211-288g17.3"
## [46] "lrmp" "syn1" "nr4a3"
## [49] "mbd3b" "atf3" "zic1"
## [52] "rab2a" "prdm8b" "etv1"
## [55] "zic2b" "foxp4" "bhlhe40"
## [58] "emx2" "zic2a" "sox3"
## [61] "kif1b" "foxp1a" "sox5"
## [64] "mef2cb" "ing4" "insm1b"
## [67] "her12" "neurod6b" "zfhx4"
## [70] "tfap2e" "jund" "hmgb1b"
## [73] "smarce1" "emx1" "etv5b"
## [76] "fezf1" "tox" "CABZ01070258.1"
## [79] "her6" "irx5a" "phf20b"
## [82] "pou3f1" "foxo1a" "neurod2"
## [85] "smyd2a" "lhx1a" "zfand5a"
## [88] "chd7" "fezf2" "lhx2b"
## [91] "zic5" "hes6" "id4"
## [94] "klf12a" "hmgb3b" "tsc22d1"
## [97] "pknox2" "her13" "erf"
## [100] "nat8l" "ZBTB7C" "znf536"
## [103] "klf2b" "max" "znf296"
## [106] "myca" "ascl1a" "ivns1abpa"
## [109] "foxg1a" "mycn" "tefa"
## [112] "olig4" "tgif1" "sim1b"
## [115] "nfil3" "onecut1" "mid1ip1b"
## [118] "pax6a" "tiparp" "scrt2"
## [121] "atf6" "pou3f3a" "tbr1a"
## [124] "her9" "rab7" "ssbp3b"
## [127] "rab11ba" "btbd6b" "hp1bp3"
## [130] "pbx1a" "etv5a" "zbtb20"
## [133] "xbp1" "tfap2a" "her7"
## [136] "edf1" "DAB2" "creb3l3l"
## [139] "mafba" "foxa1" "scrt1a"
## [142] "rcor2" "znf395b" "btbd6a"
## [145] "bcl11aa" "mxi1" "id2b"
## [148] "sim1a" "sall1a" "zdhhc2"
## [151] "mbnl2" "dbx1a" "tsc22d3"
## [154] "foxp2" "maff" "sox19a"
## [157] "rab11bb" "zfr" "foxd2"
## [160] "uncx4.1" "setd8a" "tal1"
## [163] "kdm6bb" "polr3k" "dpf1"
## [166] "skia" "srebf2" "rab11a"
## [169] "phactr2" "ssbp4" "znf346"
## [172] "roraa" "rab18a" "hif1ab"
## [175] "jph1b" "naa10" "setd3"
## [178] "lef1" "plxna4" "klhl17"
## [181] "tcf12" "zgc:171482" "znf395a"
## [184] "ccdc124" "mych" "rorab"
## [187] "dachc" "rbfox2" "zeb2a"
## [190] "raraa" "znrf1" "runx1t1"
## [193] "inab" "wdr5" "nr3c1"
## [196] "map7d1b" "tfdp2" "sox6"
## [199] "foxg1b" "tle2" "vax1"
## [202] "dacha" "mta2" "hsf2"
## [205] "nfe2" "zbtb16a" "mkxa"
## [208] "sap30bp" "hells" "irx3a"
## [211] "six3a" "mkrn1" "zbtb4"
## [214] "mrpl28" "ewsr1b" "nr2f6b"
## [217] "tcea1" "csdc2a" "kdm2aa"
## [220] "smad3a" "yy1b" "desma"
## [223] "arid3b" "nfyba" "traf4a"
## [226] "cntn1b" "slmapa" "twist2"
## [229] "nfya" "myt1la" "coq9"
## [232] "mtdhb" "spty2d1" "foxo3b"
## [235] "nr3c2" "herpud1" "plxnb1a"
## [238] "dido1" "phactr3a" "hif1al"
## [241] "nr1d2a" "zmynd8" "foxk2"
## [244] "zgc:91944" "zdhhc5a" "kdm5c"
## [247] "skib" "aebp2" "btaf1"
## [250] "mllt3" "sox1b" "hira"
## [253] "smad2" "yeats4" "ing2"
## [256] "maf1" "hopx" "IKZF2"
## [259] "zmat2" "MAF1" "zdhhc15b"
## [262] "msgn1" "nr1d2b" "six3b"
## [265] "churc1" "e2f7" "klf13"
## [268] "zic6" "id3" "klhl21"
## [271] "znf385a" "aanat2" "zgc:77880"
## [274] "ikzf5" "ncoa3" "trit1"
## [277] "cbx4" "ktn1" "supt5h"
## [280] "RHOBTB1" "smad1" "gtf2b"
## [283] "hmga1b" "klhl20" "hmg20a"
## [286] "tcf7l1b" "znhit1" "zgc:171673"
## [289] "esr1" "ubtf" "cebpg"
## [292] "cebpd" "mycb" "kdm5bb"
## [295] "atf1" "MYO5B" "zgc:65851"
## [298] "kdm6al" "rfx2" "relb"
## [301] "waslb" "stat1b" "sox19b"
## [304] "farsb" "gtf2f1" "zc3h13"
## [307] "si:dkeyp-113d7.1" "otx1a" "zdhhc6"
## [310] "pole3" "taf11" "sox2"
## [313] "smarca5" "smad7" "btbd9"
## [316] "nr1i2" "znf207b" "asf1bb"
## [319] "ssbp3a" "vox" "kdm5ba"
## [322] "keap1a" "ppardb" "stat1a"
## [325] "prdm8" "thrb" "rcc1"
## [328] "si:dkey-253d23.5" "whsc1l1" "atf7b"
## [331] "mxd1" "epas1b" "rybpa"
## [334] "ZNF608" "morc2" "dachb"
## [337] "brf1a" "znf319" "grik5"
## [340] "spopla" "mid1ip1a" "yy1a"
## [343] "rarga" "dlx6a" "hif1al2"
## [346] "zgc:114130" "ZBTB39" "dmrta2"
## [349] "zmynd11" "sp4" "zgc:66474"
## [352] "tsc22d2" "brms1" "optn"
## [355] "rfx1a" "pou2f1b" "klf11a"
## [358] "ssrp1b" "znf503" "pou6f1"
## [361] "trim33" "zc3h12b" "znf598"
## [364] "ep300a" "znf148" "gabpa"
## [367] "stat2" "pknox1.1" "NACC2"
## [370] "tfap2c" "nkx2.5" "zgc:162972"
## [373] "ash1l" "nr1d4b" "BX890616.1"
## [376] "kin" "elf2b" "atoh8"
## [379] "ing3" "plxnb3" "nfyc"
## [382] "nsd1a" "pias4a" "znf451"
## [385] "klf11b" "rel" "specc1la"
## [388] "hinfp" "pola1" "zgc:173517"
## [391] "brms1lb" "NPAS3" "ar"
## [394] "znf207a" "jarid2a" "dpf3"
## [397] "hic1" "rxraa" "znf280d"
## [400] "atoh7" "erg" "klhl32"
## [403] "rorca" "zfp36l2" "klf5b"
## [406] "ubtfl" "znf526" "foxj1a"
## [409] "nr2f6a" "mier1b" "kcmf1"
## [412] "pogzb" "noc4l" "zbtb25"
## [415] "zgc:161969" "foxn4" "mynn"
## [418] "smarcc1a" "l3mbtl1b" "sirt3"
## [421] "hmg20b" "suds3" "rnf2"
## [424] "ssrp1a" "znf142" "mxd4"
## [427] "chst11" "foxp3a" "prdm2a"
## [430] "zbtb2b" "brf1b" "taf3"
## [433] "plxnb2a" "tcea3"#plot top genes that are in AGETAZ database and vary over pseudotime
plot_cells(cds_test,
genes=deg_ids_insitu[1:36],
show_trajectory_graph = FALSE,
label_cell_groups = FALSE,
label_leaves = FALSE)
3.7 Find modules in genes that vary over time and make heatmap
#preprocess and specify number of PCs; not sure what this step does but it prevents a future error
cds_test <- preprocess_cds(cds_test, num_dim = 100)
#0e+00 1e-06 1e-05 1e-04 1e-03 1e-02 1e-01
#find gene modules
gene_modules <- find_gene_modules(cds_test[deg_ids,],
resolution=5e-03, random_seed = 1234)
#get cluster names
cluster_groups <- data.frame(cell = row.names(colData(cds_test)),
cluster_group = cds_test@clusters$UMAP[[2]])
#aggregate by modules and clusters
agg_mat2 <- aggregate_gene_expression(cds_test, gene_modules, cluster_groups)
# how many modules and clusters?
dim(agg_mat2)## [1] 37 9#make rownames of agg_mat2 the names of the modules
row.names(agg_mat2) <- paste0("Module ", row.names(agg_mat2))
#make heatmap with modules and clusters
pheatmap::pheatmap(agg_mat2,
scale="column",
treeheight_row = 0,
treeheight_col = 0,
clustering_method="ward.D2")
#export gene modules
write.csv(gene_modules, "gene_modules.csv")This resolution produces 37 modules.
Progenitors: 10, 20, 35 Pallium 02/03: 11?, 19, 23, 25 Pallium 05/06/07: 9,13,14,16,17
3.8 Plot expression of genes in modules on UMAP
#plot all modules
plot_cells(cds_test,
genes=gene_modules,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
3.9 Look at top genes of select progenitor modules
#plot some progenitor modules
plot_cells(cds_test,
genes=gene_modules %>% filter(module %in% c(35, 20, 10)),
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#what are the top genes marking these modules?
module35genes <- unlist((gene_modules %>% filter(module %in% c(35)))[1:16,1])
#plot top 16 genes of module 35
plot_cells(cds_test,
genes=module35genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 20 genes
module20genes <- unlist((gene_modules %>% filter(module %in% c(20)))[1:16,1])
#plot top 16 genes of module 20
plot_cells(cds_test,
genes=module20genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 10 genes
module10genes <- unlist((gene_modules %>% filter(module %in% c(10)))[1:16,1])
#plot top 16 genes of module 20
plot_cells(cds_test,
genes=module10genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
These look okay, but module 10 is all ribosomal genes? I guess these cells are making proteins?
3.10 Look at top genes of select Pallium 02/03 modules
#plot some pallium modules
plot_cells(cds_test,
genes=gene_modules %>% filter(module %in% c(11, 19, 23, 25)),
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 11 genes
module11genes <- unlist((gene_modules %>% filter(module %in% c(11)))[1:16,1])
#plot top 16 genes of module 11
plot_cells(cds_test,
genes=module11genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 19 genes
module19genes <- unlist((gene_modules %>% filter(module %in% c(19)))[1:16,1])
#plot top 16 genes of module 19
plot_cells(cds_test,
genes=module19genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 23 genes
module23genes <- unlist((gene_modules %>% filter(module %in% c(23)))[1:16,1])
#plot top 16 genes of module 23
plot_cells(cds_test,
genes=module23genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 25 genes
module25genes <- unlist((gene_modules %>% filter(module %in% c(25)))[1:16,1])
#plot top 16 genes of module 25
plot_cells(cds_test,
genes=module25genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
Module 11 isn’t great. Module 19 is pretty good. Should check expression of a few module 19 genes (cbln1, rspo2) in the larger dataset. Some of the genes in module 23 look pretty specific but mostly low expression. Module 25 is also decent and finds pbx3b, lhx5, and lhx9. Overall 19 and 25 look the best for pallium02/03.
3.11 Look at top genes of select Pallium 05/06/07 modules
#plot some pallium modules
plot_cells(cds_test,
genes=gene_modules %>% filter(module %in% c(9,13,14,16,17)),
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 9 genes
module9genes <- unlist((gene_modules %>% filter(module %in% c(9)))[1:16,1])
#plot top 16 genes of module 9
plot_cells(cds_test,
genes=module9genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 13 genes
module13genes <- unlist((gene_modules %>% filter(module %in% c(13)))[1:16,1])
#plot top 16 genes of module 13
plot_cells(cds_test,
genes=module13genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 14 genes
module14genes <- unlist((gene_modules %>% filter(module %in% c(14)))[1:16,1])
#plot top 16 genes of module 14
plot_cells(cds_test,
genes=module14genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 16 genes
module16genes <- unlist((gene_modules %>% filter(module %in% c(16)))[1:16,1])
#plot top 16 genes of module 16
plot_cells(cds_test,
genes=module16genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
#module 17 genes
module17genes <- unlist((gene_modules %>% filter(module %in% c(17)))[1:16,1])
#plot top 16 genes of module 17
plot_cells(cds_test,
genes=module17genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
None of these modules really have genes expressed in all of pallium 05 and 06. Module 9 might be the best. Might want to look up bhlhe23(module 16), nell2b (module 9), ets2 (module 9), fam43b(module 9), si:ch211-51a6.2 (module 9).
4 Back to Seurat
4.1 Add pseudotime information to original Seurat object & plot
#load original object
load(file = "forebrain_all_final.RObj")
#add pseudotime to seurat object
cds_test$monocle3_pseudotime <- pseudotime(cds_test)
data.pseudo <- as.data.frame(colData(cds_test))
forebrain.integrated$pseudotime <- pseudotime(cds_test)
#plot pseudotime on seurat object
FeaturePlot(forebrain.integrated, features = "pseudotime")
4.2 Make plots with some possibly interesting genes
#module 35
FeaturePlot(forebrain.integrated, features = module35genes)
#module 20
FeaturePlot(forebrain.integrated, features = module20genes)
#module 19
FeaturePlot(forebrain.integrated, features = module19genes)
#module 25
FeaturePlot(forebrain.integrated, features = module25genes)
#module 9
FeaturePlot(forebrain.integrated, features = module9genes)
4.3 Make ridge plots
Could make these for specific genes/clusters if we want to. This is just an example with genes from modules 35, 20, 19, and 25.
RidgePlot(forebrain.integrated, features = c('cbln1'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
RidgePlot(forebrain.integrated, features = c('her4.2'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
RidgePlot(forebrain.integrated, features = c('neurod4'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
RidgePlot(forebrain.integrated, features = c('lhx9'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
RidgePlot(forebrain.integrated, features = c('pbx3b'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
RidgePlot(forebrain.integrated, features = c('elavl4'), sort = F, idents = c('Committed_pallium_precursors', 'Pallium_IN02', 'Pallium_02', "Pallium_03"))
4.4 Testing GO enrichment on modules
library(gprofiler2) #tools for accessing the GO enrichment results using g:Profiler web resources
#testing with ribosomal module
module10genes <- unlist((gene_modules %>% filter(module %in% c(10)))[,1])
# use the 'gost' function from the gprofiler2 package to run GO enrichment analysis
gost.res <- gost(module10genes, organism = "drerio", correction_method = "fdr")
# produce an interactive manhattan plot of enriched GO terms
gostplot(gost.res, interactive = T, capped = T) #set interactive=FALSE to get plot for publications#module 25
modulegenes <- unlist((gene_modules %>% filter(module %in% c(25)))[,1])
gost.res <- gost(modulegenes, organism = "drerio", correction_method = "fdr")
gostplot(gost.res, interactive = T, capped = T) #module 9
modulegenes <- unlist((gene_modules %>% filter(module %in% c(9)))[,1])
gost.res <- gost(modulegenes, organism = "drerio", correction_method = "fdr")
gostplot(gost.res, interactive = T, capped = T)#module 19
modulegenes <- unlist((gene_modules %>% filter(module %in% c(19)))[,1])
gost.res <- gost(modulegenes, organism = "drerio", correction_method = "fdr")
gostplot(gost.res, interactive = T, capped = T)4.5 Session info
Packages and versions necessary to reproduce the results in this report.
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux
##
## Matrix products: default
## BLAS/LAPACK: /data/rc/apps/rc/software/FlexiBLAS/3.0.4-GCC-10.3.0/lib64/libflexiblas.so.3.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gprofiler2_0.2.1 patchwork_1.1.2
## [3] SeuratWrappers_0.3.1 monocle3_1.3.1
## [5] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
## [7] GenomicRanges_1.50.2 GenomeInfoDb_1.34.6
## [9] IRanges_2.32.0 S4Vectors_0.36.1
## [11] MatrixGenerics_1.10.0 matrixStats_0.63.0
## [13] Biobase_2.58.0 BiocGenerics_0.44.0
## [15] SeuratObject_4.1.3 Seurat_4.1.1
## [17] forcats_0.5.2 stringr_1.5.0
## [19] dplyr_1.0.10 purrr_1.0.0
## [21] readr_2.1.3 tidyr_1.2.1
## [23] tibble_3.1.8 ggplot2_3.4.0
## [25] tidyverse_1.3.2 knitr_1.41
## [27] tinytex_0.43 rmarkdown_2.19
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.26
## [3] R.utils_2.12.2 tidyselect_1.2.0
## [5] lme4_1.1-31 htmlwidgets_1.6.1
## [7] grid_4.2.0 Rtsne_0.16
## [9] munsell_0.5.0 codetools_0.2-18
## [11] ica_1.0-3 future_1.30.0
## [13] miniUI_0.1.1.1 withr_2.5.0
## [15] spatstat.random_3.1-0.008 colorspace_2.0-3
## [17] progressr_0.12.0 highr_0.10
## [19] rstudioapi_0.14 ROCR_1.0-11
## [21] tensor_1.5 listenv_0.9.0
## [23] labeling_0.4.2 GenomeInfoDbData_1.2.9
## [25] polyclip_1.10-4 pheatmap_1.0.12
## [27] farver_2.1.1 parallelly_1.33.0
## [29] vctrs_0.5.1 generics_0.1.3
## [31] xfun_0.36 timechange_0.1.1
## [33] R6_2.5.1 ggbeeswarm_0.7.1
## [35] rsvd_1.0.5 bitops_1.0-7
## [37] spatstat.utils_3.0-1 cachem_1.0.6
## [39] DelayedArray_0.24.0 assertthat_0.2.1
## [41] promises_1.2.0.1 scales_1.2.1
## [43] googlesheets4_1.0.1 beeswarm_0.4.0
## [45] gtable_0.3.1 globals_0.16.2
## [47] goftest_1.2-3 rlang_1.0.6
## [49] splines_4.2.0 lazyeval_0.2.2
## [51] gargle_1.2.1 spatstat.geom_3.0-4.001
## [53] broom_1.0.2 BiocManager_1.30.19
## [55] yaml_2.3.6 reshape2_1.4.4
## [57] abind_1.4-5 modelr_0.1.10
## [59] crosstalk_1.2.0 backports_1.4.1
## [61] httpuv_1.6.7 tools_4.2.0
## [63] ellipsis_0.3.2 spatstat.core_2.4-4.010
## [65] jquerylib_0.1.4 RColorBrewer_1.1-3
## [67] proxy_0.4-27 ggridges_0.5.4
## [69] Rcpp_1.0.9 plyr_1.8.8
## [71] sparseMatrixStats_1.10.0 zlibbioc_1.44.0
## [73] RCurl_1.98-1.9 rpart_4.1.19
## [75] deldir_1.0-6 viridis_0.6.2
## [77] pbapply_1.6-0 cowplot_1.1.1
## [79] zoo_1.8-11 grr_0.9.5
## [81] haven_2.5.1 ggrepel_0.9.2
## [83] cluster_2.1.4 fs_1.5.2
## [85] magrittr_2.0.3 data.table_1.14.6
## [87] scattermore_0.8 lmtest_0.9-40
## [89] reprex_2.0.2 RANN_2.6.1
## [91] googledrive_2.0.0 fitdistrplus_1.1-8
## [93] hms_1.1.2 mime_0.12
## [95] evaluate_0.19 xtable_1.8-4
## [97] readxl_1.4.1 gridExtra_2.3
## [99] compiler_4.2.0 KernSmooth_2.23-20
## [101] crayon_1.5.2 minqa_1.2.5
## [103] R.oo_1.25.0 htmltools_0.5.4
## [105] mgcv_1.8-41 later_1.3.0
## [107] tzdb_0.3.0 lubridate_1.9.0
## [109] DBI_1.1.3 dbplyr_2.2.1
## [111] MASS_7.3-58.1 boot_1.3-28.1
## [113] leidenbase_0.1.14 Matrix_1.5-3
## [115] cli_3.5.0 R.methodsS3_1.8.2
## [117] parallel_4.2.0 igraph_1.3.5
## [119] pkgconfig_2.0.3 sp_1.5-1
## [121] terra_1.6-47 plotly_4.10.1
## [123] spatstat.sparse_3.0-0 xml2_1.3.3
## [125] vipor_0.4.5 bslib_0.4.2
## [127] XVector_0.38.0 rvest_1.0.3
## [129] digest_0.6.31 sctransform_0.3.5
## [131] RcppAnnoy_0.0.20 spatstat.data_3.0-0
## [133] cellranger_1.1.0 leiden_0.4.3
## [135] uwot_0.1.14 DelayedMatrixStats_1.20.0
## [137] shiny_1.7.4 nloptr_2.0.3
## [139] lifecycle_1.0.3 nlme_3.1-161
## [141] jsonlite_1.8.4 viridisLite_0.4.1
## [143] fansi_1.0.3 pillar_1.8.1
## [145] lattice_0.20-45 ggrastr_1.0.1
## [147] fastmap_1.1.0 httr_1.4.4
## [149] survival_3.4-0 glue_1.6.2
## [151] remotes_2.4.2 png_0.1-8
## [153] stringi_1.7.8 sass_0.4.4
## [155] irlba_2.3.5.1 future.apply_1.10.0